library(DAseq)
library(Seurat)
Attaching SeuratObject
library(ggplot2)
library(patchwork)
library(stringr)
# set the python path
python2use <- "~/anaconda3/bin/python3"

my.samplename <- "Gg_ctrl_poly_int"

my.se <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_seurat_070223.rds")

DimPlot(my.se, reduction = "tsne", label = TRUE)

DA scores

cell.labels <- unname(substr(my.se$orig.ident, 1, 4))

resp <- "poly"
non_resp <- "ctrl"

# Get the DA scores for all cells
da.cells <- getDAcells(
  X = my.se@reductions$pca@cell.embeddings[, 1:30],
  cell.labels = cell.labels,
  labels.1 = resp, #responder labels (blue)
  labels.2 = non_resp, # nonresponder labels (red)
  k.vector = seq(50, 500, 50),
  plot.embedding = my.se@reductions$tsne@cell.embeddings
)
Calculating DA score vector.
Running GLM.
Test on random labels.
Setting thresholds based on permutation
# same for UMAP
da.cells.u <- getDAcells(
  X = my.se@reductions$pca@cell.embeddings[, 1:30],
  cell.labels = cell.labels,
  labels.1 = resp, #responder labels (blue)
  labels.2 = non_resp, # nonresponder labels (red)
  k.vector = seq(50, 500, 50),
  plot.embedding = my.se@reductions$umap@cell.embeddings
)
Calculating DA score vector.
Running GLM.
Test on random labels.
Setting thresholds based on permutation

With a random permutation ofn the labels, a null distribution is generated. Min and max of the permutation DA set the default threshold.

# look at output
str(da.cells[1:2])
List of 2
 $ cell.idx: int [1:15537] 1 2 3 4 5 6 7 8 9 10 ...
 $ da.ratio: num [1:15537, 1:10] -0.4171 0.0939 -0.0272 0.0939 0.0939 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:10] "50" "100" "150" "200" ...
# plot the vis of the DA scores
p1 <- da.cells$pred.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))

p2 <- DimPlot(my.se,
        ncol = 4,
        reduction = "tsne",
        pt.size = 0.1,
        split.by = "orig.ident")

p3 <- DimPlot(my.se, 
        group.by = "orig.ident",
        reduction = "tsne",
        pt.size = 1)

p4 <- da.cells.u$pred.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))

p5 <- DimPlot(my.se,
        ncol = 4,
        reduction = "umap",
        pt.size = 0.1,
        split.by = "orig.ident")

p6 <- DimPlot(my.se, 
        group.by = "orig.ident",
        reduction = "umap",
        pt.size = 1)

p2 

p3 + p1

p5

p6 + p4

# plot randomisation of DA scores
da.cells$rand.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue/down) and Non-responder \"", non_resp,"\" (red/up)"))

With updateDAcells we can apply our own DA threshold to be more (or less) stringent.

# update DA cells with custom threshold
da.cells <- updateDAcells(
  X = da.cells,
  pred.thres = c(-0.7,0.7),
  plot.embedding = my.se@reductions$tsne@cell.embeddings
)

# check output
da.cells$da.cells.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))

da.cells.u <- updateDAcells(
  X = da.cells,
  pred.thres = c(-0.7,0.7),
  plot.embedding = my.se@reductions$umap@cell.embeddings
)

# check output
da.cells.u$da.cells.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))

Regions

# get the contiguous regions of cells
da.regions <- getDAregion(
  X = my.se@reductions$pca@cell.embeddings[, 1:20],
  da.cells = da.cells,
  cell.labels = cell.labels,
  labels.1 = resp,
  labels.2 = non_resp,
  resolution = 0.01,
  plot.embedding = my.se@reductions$tsne@cell.embeddings,
)
Using min.cell = 50
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: k.param set larger than number of cells. Setting k.param to number of
cells - 1.
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Removing 4 DA regions with cells < 50.
str(da.regions[1:2])
List of 2
 $ cell.idx       : int [1:15537] 1 2 3 4 5 6 7 8 9 10 ...
 $ da.region.label: num [1:15537] 0 0 0 0 0 0 0 0 0 0 ...
da.regions$da.region.plot

# get the contiguous regions of cells (UMAP)
da.regions.u <- getDAregion(
  X = my.se@reductions$pca@cell.embeddings[, 1:20],
  da.cells = da.cells,
  cell.labels = cell.labels,
  labels.1 = resp,
  labels.2 = non_resp,
  resolution = 0.01,
  plot.embedding = my.se@reductions$umap@cell.embeddings,
)
Using min.cell = 50
Warning: The following arguments are not used: row.names

Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: k.param set larger than number of cells. Setting k.param to number of
cells - 1.
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Removing 4 DA regions with cells < 50.
da.regions.u$da.region.plot

# add da regions as meta data to plot them
my.se$da.regions <- da.regions$da.region.label

Export

# cell x gene matrix
expr <- as.matrix(my.se[["RNA"]]@data)

STG.markers <- STGmarkerFinder(
  X = expr,
  da.regions = da.regions,
  lambda = 1.5, n.runs = 5, return.model = T,
  python.use = python2use
)

# wilcoxon rank sum (MAST does not work)
Seurat.markers <-SeuratMarkerFinder(
  object = my.se,
  features = my.se[["integrated"]]@var.features,
  da.slot = "da.regions",
  da.regions.to.run = 1:(length(table(da.regions$da.region.label))-1),
  min.pct = 0.25,
  logfc.threshold =  0.5,
  return.thresh = 0.05,
  assay = "RNA"
  )

Save the output:

saveRDS(da.cells, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_cells.rds"))
saveRDS(da.regions, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_regions.rds"))
saveRDS(da.cells.u, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_cells_umap.rds"))
saveRDS(da.regions.u, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_regions_umap.rds"))
saveRDS(STG.markers, paste0("~/spinal_cord_paper/output/",my.samplename,"_STG_markers.rds"))
saveRDS(Seurat.markers, paste0("~/spinal_cord_paper/output/",my.samplename,"_Seurat_markers.rds"))
pdf(paste0("~/spinal_cord_paper/figures/",my.samplename,"_DAseq.pdf"), width = 10, height = 10)
DimPlot(
  my.se,
  group.by = "da.regions",
  reduction = "tsne",
  cols = c(
    "grey",
    rainbow(
      length(table(da.regions$da.region.label))-1
      )
    ),
  label = TRUE,
  repel = TRUE,
  label.size = 6
  )+
  ggtitle(paste0(my.samplename," DA regions"))+
  theme(plot.title = element_text(hjust = 0.5))
DimPlot(
  my.se,
  group.by = "da.regions",
  reduction = "umap",
  cols = c(
    "grey",
    rainbow(
      length(table(da.regions$da.region.label))-1
      )
    ),
  label = TRUE,
  repel = TRUE,
  label.size = 6
  )+
  ggtitle(paste0(my.samplename," DA regions"))+
  theme(plot.title = element_text(hjust = 0.5))
p2
p1
p5
p4
da.cells$rand.plot
da.cells$da.cells.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
da.cells.u$da.cells.plot + 
  ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
dev.off()
png 
  2 
# Date and time of Rendering
Sys.time()
[1] "2023-02-27 12:12:45 CET"
sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /scicore/soft/apps/OpenBLAS/0.3.1-GCC-7.3.0-2.30/lib/libopenblas_sandybridgep-r0.3.1.so

locale:
[1] en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stringr_1.4.0      patchwork_1.1.1    ggplot2_3.3.3      SeuratObject_4.0.2
[5] Seurat_4.0.5       DAseq_1.0.0       

loaded via a namespace (and not attached):
  [1] plyr_1.8.6            igraph_1.2.6          lazyeval_0.2.2       
  [4] sp_1.4-5              splines_4.1.0         listenv_0.8.0        
  [7] scattermore_0.7       digest_0.6.27         foreach_1.5.1        
 [10] htmltools_0.5.1.1     fansi_0.5.0           magrittr_2.0.1       
 [13] tensor_1.5            cluster_2.1.2         ROCR_1.0-11          
 [16] limma_3.48.0          recipes_0.1.16        globals_0.16.2       
 [19] gower_0.2.2           matrixStats_0.58.0    spatstat.sparse_3.0-0
 [22] colorspace_2.0-1      ggrepel_0.9.1         xfun_0.34            
 [25] dplyr_1.0.10          jsonlite_1.7.2        spatstat.data_3.0-0  
 [28] survival_3.2-11       zoo_1.8-9             iterators_1.0.13     
 [31] glue_1.6.2            polyclip_1.10-0       gtable_0.3.0         
 [34] ipred_0.9-11          leiden_0.3.9          future.apply_1.7.0   
 [37] shape_1.4.6           abind_1.4-5           scales_1.1.1         
 [40] DBI_1.1.1             miniUI_0.1.1.1        Rcpp_1.0.7           
 [43] viridisLite_0.4.0     xtable_1.8-4          reticulate_1.22      
 [46] spatstat.core_2.1-2   proxy_0.4-25          stats4_4.1.0         
 [49] lava_1.6.9            prodlim_2019.11.13    glmnet_4.1-1         
 [52] htmlwidgets_1.5.3     httr_1.4.2            RColorBrewer_1.1-2   
 [55] ellipsis_0.3.2        ica_1.0-2             farver_2.1.0         
 [58] pkgconfig_2.0.3       nnet_7.3-16           sass_0.4.0           
 [61] uwot_0.1.10           deldir_1.0-6          here_1.0.1           
 [64] utf8_1.2.1            caret_6.0-88          tidyselect_1.2.0     
 [67] labeling_0.4.2        rlang_1.0.6           reshape2_1.4.4       
 [70] later_1.2.0           munsell_0.5.0         tools_4.1.0          
 [73] cli_3.4.1             generics_0.1.3        ggridges_0.5.3       
 [76] evaluate_0.20         fastmap_1.1.0         yaml_2.2.1           
 [79] goftest_1.2-2         ModelMetrics_1.2.2.2  knitr_1.41           
 [82] fitdistrplus_1.1-6    purrr_0.3.4           RANN_2.6.1           
 [85] pbapply_1.4-3         future_1.30.0         nlme_3.1-152         
 [88] mime_0.10             compiler_4.1.0        plotly_4.10.0        
 [91] png_0.1-7             e1071_1.7-7           spatstat.utils_3.0-1 
 [94] tibble_3.1.8          bslib_0.2.5.1         stringi_1.6.2        
 [97] highr_0.9             lattice_0.20-44       Matrix_1.3-3         
[100] vctrs_0.5.1           pillar_1.8.1          lifecycle_1.0.3      
[103] spatstat.geom_3.0-3   lmtest_0.9-38         jquerylib_0.1.4      
[106] RcppAnnoy_0.0.19      data.table_1.14.0     cowplot_1.1.1        
[109] irlba_2.3.3           httpuv_1.6.1          R6_2.5.0             
[112] promises_1.2.0.1      KernSmooth_2.23-20    gridExtra_2.3        
[115] parallelly_1.33.0     codetools_0.2-18      MASS_7.3-54          
[118] assertthat_0.2.1      rprojroot_2.0.2       withr_2.4.2          
[121] sctransform_0.3.3     mgcv_1.8-35           parallel_4.1.0       
[124] grid_4.1.0            rpart_4.1-15          timeDate_3043.102    
[127] tidyr_1.1.3           class_7.3-19          rmarkdown_2.17       
[130] Rtsne_0.15            pROC_1.17.0.1         shiny_1.6.0          
[133] lubridate_1.7.10